Experimental enforcement of sex-limited autosome inheritance does not reveal intralocus sexual conflict

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = FALSE)

Load packages

library(tidyverse) # tidy style coding
library(png) # to load images
library(grid) # to plot images
library(brms) # Bayesian models
library(loo) # to use information criteria in brms models
library(dagitty) # to build DAGs
library(ggdag) # to tidy DAGs
library(tidybayes) # Bayesian aesthetics
library(MetBrewer) # colours
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(pander) # for simpler tables
library(stickylabeller) # for labelling facets in ggplot
library(ggtext) # for markdown notation in labs
library(geomtextpath) # for labelling contours

Supplementary methods

\(~\)

LH_M culturing

We maintained this LH_M stock in our laboratory for 32 generations prior to creating the genotypes required for experimental evolution. We cultured LH_M at 25C, with a 16-8 light-dark cycle and reared in vials (95mm x 25mm) on a corn-meal, yeast and dextrose-based diet (recipe in Table S1; ~8cm3 of food medium per vial) supplemented with dried yeast, at a population size of at least 800 breeding individuals across 25 vials (16 flies of each sex per vial, following Rice et al. 2005). Each generation begins by pooling the offspring produced across the 25 vials and randomly assorting 16 female-male pairs to 25 new vials. We then allow these breeding individuals 48 hours to interact and mate, before transferring them to another set of new vials. After 24 hours of egg-laying, we discard all adults, and allow juveniles 12 days to compete for resources, pupate and eclose as adults. We then iteratively repeat this process each generation to maintain the population.

\(~\)

Creating the genotypes used for experimental evolution

img <- readPNG("Figure_S1.png")
 
grid.raster(img)

Figure S1. Crossing scheme used to integrate the GFP constructs and ap^XA marked translocated second and third chromosome balancers into the LH_M genetic background. We replicated the crosses 12 times to supply the flies used in generation zero of experimental evolution; 6 times using the Ubi GFP construct and 6 times with the 3xP GFP construct. We performed each cross across 25 vials, with 16 females and 16 males in each, to preserve genetic variation in our evolving populations. G = generation.

\(~\)

Table S1. Recipe for food medium used in our experiment. The provided quantities make ~ 1 litre of food.

tibble("Ingredients" = c("Soy flour", "Cornmeal", "Yeast", "Dextrose", "Agar", "Water", "Tegosept", "Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL)"),
       "Quantity" = c("20 g", "73 g", "35 g", "75 g", "6 g", "1000 mL", "17 mL", "14 mL")) %>% 
  pander(split.cell = 40, split.table = Inf)

Ingredients	Quantity
Soy flour	20 g
Cornmeal	73 g
Yeast	35 g
Dextrose	75 g
Agar	6 g
Water	1000 mL
Tegosept	17 mL
Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL)	14 mL

\(~\)

The recombination compartment

\(~\)

The female recombination compartment

When we initiated the first generation of experimental evolution for each of the sex-limited populations, we placed 12 females with the FLA/FLA genotype into individual food vials, each containing a FLA/ap^XA male. We allowed them ~24 hours to mate, then discarded the males and pooled the 12 females into a single vial. After ~72 hours we discarded the females and allowed their offspring to develop. The aim of this breeding design was to minimise selection acting on the FLA/ap^XA males, since each male simply needed to survive and then fertilise one randomly assigned virgin female. From the progeny, we collected 24 virgin female offspring with the genotype FLA/FLA, where recombination had occurred between the homologous FLA chromosomes. Half of these females were used to establish the next iteration of the recombination compartment, where they were again individually mated to 12 FLA/ap^XA males sourced from the progeny of the main population. The 12 remaining FLA/FLA females, which carried one set of recombined FLA autosomes inherited from their mother, were randomly substituted for 12 females in the main breeding population. Therefore, 6% of the females in each generation came from the recombination compartment (Figure 1).

The male recombination compartment

At the start of the first generation for each sex-limited population, we placed 12 females with genotype MLA/ap^XA into a single food vial along with 12 MLA/ap^XA males. We allowed ~24 hours for the males to compete for fertilisations, after which we discarded the males and moved the females into individual food vials to oviposit. We collected one female and one male MLA/MLA offspring from each of the oviposition vials. This rearing protocol minimised selection on females, since females were reared mostly in a competition- and harassment-free environment, and we equalised fitness by collecting a standard number of progeny from each female. The 12 males were randomly substituted for 12 males in the main breeding population breeding population (carrying recombined MLA autosomes inherited from their mother), and the 12 females were used in the next iteration of the male recombination compartment. Therefore, 6% of the males in each generation came from the recombination compartment (Figure 1).

We also included a female and male recombination compartment for each of the four control populations. The process was identical as for the sex-limited populations, except that control autosomes can be substituted in the above descriptions for female- or male-limited autosomes.

\(~\)

Modelling the effect of the relaxed selection phase

\(~\)

Relaxation of sex-limited inheritance may have allowed sexually antagonistic allele frequencies to rebound back to their original pre-treatment equilibria, nullifying any effect of our experiment. To assess the effect of the relaxed selection phase on our experiment, we built a single-locus, two-allele population genetic model. The model contains three phases: in phase one, we recapitulated prior theory to find the conditions where intralocus sexual conflict is expected to maintain polymorphism. These conditions provide a baseline expectation for the characteristics of the alleles targeted by our expeirmental evolution treatments. Then, in phase two, we modelled how allele frequencies are expected to change under 20 generations of sex-limited evolution. Finally, in phase three, we modelled rebound change in allele frequency when sex-limited evolution was relaxed.

In each phase, the model follows the dynamics of a single autosomal locus with two alleles, \(A_f\) and \(A_m\), which have frequencies p and 1 - p, respectively. For simplicity, we assume an infinite population size, discrete, non-overlapping generations and panmixia with respect to the focal locus. Following previous models of sexual antagonism (Kidwell et al, 1977; Rice, 1984; Connallon et al, 2009; Patten, et al, 2010), we create intralocus sexual conflict by coding the \(A_f\) allele to have greater fitness than the \(A_m\) allele when expressed in females, while the \(A_m\) allele has greater fitness than the \(A_f\) allele when expressed in males. The fitness benefits of being homozygous for the ‘correct’ allele in the ‘correct’ sex are represented by \(s_f\) and \(s_m\) (previous models have generally considered the cost of being homozygous for the wrong allele; we make this small change to improve the conceptual flow between the three phases of our model). In heterozygotes, these benefits are moderated by the sex-specific dominance coefficients, \(h_f\) and \(h_m\) (following Curtsinger et al 1994). When \(h_f = 0\), the \(A_f\) allele is fully recessive when expressed in females, and fully dominant when \(h_f = 1\). Modelling sex-specific dominance allows us to explore dominance reversals, where the \(A_f\) allele is dominant in one sex but recessive in the other, a mechanism that has previously been shown to facilitate balancing selection (reviewed in Connallon and Chenoweth, 2019). Table S2 shows the sex-specific fitness of each genotype.

tibble(Sex = c("Female", "Male"),
       "$A_fA_f$" = c("$1 + s_f$", "1"),
       "$A_fA_m$" = c("$1 + s_fh_f$", "$1 + s_mh_m$"),
       "$A_mA_m$" = c("1", "$1 + s_m$")) %>% 
kbl(format = "html", booktabs = T, escape = FALSE,
    caption = "<b>Table S2<b>. Effects of genotype on sex-specific fitness") %>%
  kable_styling() %>% 
  add_header_above(c(" " = 1, "Genotypes" = 3))

Table S2. Effects of genotype on sex-specific fitness

	Genotypes
Sex	\(A_fA_f\)	\(A_fA_m\)	\(A_mA_m\)
Female	\(1 + s_f\)	\(1 + s_fh_f\)	1
Male	1	\(1 + s_mh_m\)	\(1 + s_m\)

\(~\)

Phase 1: modelling polymorphism under intralocus sexual conflict

Calculating the per generation change in allele frequency

To find regions of stable polymorphism created by intralocus sexual conflict, we first find the within generation change in allele frequency following selection. For females the recursion equation is

\[\begin{equation} p'_f = \frac{p^2(1 + s_f) + \frac{1}{2}[2p(1-p)] (1 + s_fh_f)}{p^2(1 + s_f) + [2p(1-p)] (1 + s_fh_f) + (1 - p)^2} \end{equation}\]

and for males

\[\begin{equation} p'_m = \frac{p^2 + \frac{1}{2}[2p(1-p)] (1 + s_m h_m)}{p^2 + [2p(1-p)] (1 + s_mh_m) + (1 - p)^2(1 + s_m)} \end{equation}\]

Assuming an even primary sex ratio and mendelian inheritance, the \(A_f\) allele frequency in the following generation is

\[\begin{equation} p = \frac{p'_f + p'_m}{2} \end{equation}\]

We then used equation 3 to numerically determine the selective conditions required for polymorphism at this locus. We explored a parameter space of varying values of s_f, s_m, h_f and h_m, as these parameters have previously been found to be important factors in determining polymorphism (Kidwell et al, 1977; Curtsinger et al, 1994). We set the initial frequency of the \(A_f\) allele to 0.01 and iterated until 10000 generations had elapsed.

Test that the equation is correct by inputing some values for each parameter:

s_m <- 0.7

s_f <- 0.6

p <- 0.01

h_m <- 0.5

h_f <- 0.5

# female freq after selection

f_2 <- ((p^2)*(1 + s_f) + 0.5*(2*p*(1 - p))*(1 + s_f*h_f))/(
 (p^2)*(1 + s_f) + (2*p*(1 - p))*(1 + s_f*h_f) + (1-p)^2
)

# male freq after selection

m_2 <- ((p^2) + 0.5*(2*p*(1 - p))*(1 + s_m*h_m))/(
 (p^2) + (2*p*(1 - p))*(1 + s_m*h_m) + ((1-p)^2)*(1 + s_m)
)

(f_2 + m_2)/2

## [1] 0.01045281

Build the diploid simulation function

Now that we can calculate the A_f allele frequency across generations, we can find conditions where stable polymorphism is possible. To do this we build a simulation where the A_f allele is initially rare, but changes in frequency over the course of 10000 generations.

run_diploid_simulation <- function(gens_1, row, parameters){
  
  
  generations_1 <- gens_1
  s_f <- parameters$s_f[row]
  s_m <- parameters$s_m[row]
  h_f <- parameters$h_f[row]
  h_m <- parameters$h_m[row]
  initial_freq_F <- parameters$initial_freq_F[row]
  
  
  # Make the initial population
  pop <- 
    tibble(prop_1 = c(initial_freq_F),
           generations_1 = gens_1)
  
  
  current_generation <- 1
  
  # Iterate over generations
  
  while(current_generation <= generations_1){
    pop <- pop %>%
      mutate(# here's the effect of selection
        f_freq_after_selection = ((prop_1^2)*(1 + s_f) + 0.5*(2*prop_1*(1 - prop_1))*(1 + s_f*h_f))/(
          (prop_1^2)*(1 + s_f) + (2*prop_1*(1 - prop_1))*(1 + s_f*h_f) + (1-prop_1)^2),
        
        m_freq_after_selection = ((prop_1^2) + 0.5*(2*prop_1*(1 - prop_1))*(1 + s_m*h_m))/(
          (prop_1^2) + (2*prop_1*(1 - prop_1))*(1 + s_m*h_m) + ((1-prop_1)^2)*(1 + s_m)),
        
        prop_1 = (f_freq_after_selection + m_freq_after_selection)/2,
        generations_1 = gens_1)
    
    # Otherwise, increment generation counter by 1 and return to top of the 'while' loop
    
    current_generation <- current_generation + 1
  }
  
  pop %>% bind_cols(parameters[row, ]) # returns the final pop
  
}

Define the relevant parameter space

We vary \(s_f\), \(s_m\), \(h_f\) and \(h_m\) as these have previously been found to be important factors in determining polymorphism.

resolution <- 21

equilibrium_parameters <- 
  expand_grid(
    s_f = c(seq(0, 1, length = resolution*4)),
    s_m = c(seq(0, 1, length = resolution*4)),
    h_f = c(0, 0.25, 0.5, 0.75, 1),
    h_m = c(0, 0.25, 0.5, 0.75, 1),
    initial_freq_F =  0.01) %>%  
  mutate(parameter_space_ID = row_number())

Run the phase 1 simulation

if(!file.exists("polymorphism_results.csv")){
  polymorphism_results <- run_diploid_simulation(gens_1 = 10000, 1:nrow(equilibrium_parameters), equilibrium_parameters)
  write_csv(polymorphism_results, file = "polymorphism_results.csv")
} else polymorphism_results <- read_csv("polymorphism_results.csv")

Plot the results

pal1 <- met.brewer("Hiroshige", n=20, direction = -1)


polymorphism_results %>%
  mutate(Dominance = case_when(h_f == 0.5 & h_m == 0.5 ~ "Additive",
                               h_f == 1 & h_m == 0 ~ "Dominant",
                               h_f == 0 & h_m == 1 ~ "Recessive",
                               h_f == 0.25 & h_m == 0.25 ~ 
                                 "Deleterious dominance reversal (partial)",
                               h_f == 0.75 & h_m == 0.75 ~ 
                                 "Beneficial dominance reversal (partial)",
                               h_f == 1 & h_m == 1 ~ 
                                 "Beneficial dominance reversal (complete)",
                               .default = NA)) %>% 
  filter(Dominance == "Additive" | 
           Dominance == "Dominant" |
           Dominance == "Recessive" |
           Dominance == "Beneficial dominance reversal (partial)" |
           Dominance == "Beneficial dominance reversal (complete)" |
           Dominance == "Deleterious dominance reversal (partial)") %>% 
  mutate(Dominance = fct_relevel(Dominance, 
                                 "Additive",
                                 "Dominant",
                                 "Recessive",
                                 "Deleterious dominance reversal (partial)", 
                                 "Beneficial dominance reversal (partial)", 
                                 "Beneficial dominance reversal (complete)")) %>% 
  mutate(prop_1 = round(prop_1, 3)) %>% 
  ggplot(aes(x = s_f, y = s_m)) +
  geom_blank() +
  geom_raster(aes(fill = prop_1), alpha = 0.75) +
  scale_fill_gradientn(colours = pal1, breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), 
                       limits = c(0, 1)) +
  stat_contour(aes(z = prop_1), colour = "black", breaks = c(0.01, 0.99),
               linetype = 2, linewidth = 1) +
  #stat_contour(aes(z = prop_1), colour = "black", breaks = c(0.2, 0.4, 0.6, 0.8),
  #            linetype = 1, linewidth = 0.5) +
  geom_textcontour(aes(z = prop_1), breaks = c(0.2, 0.4, 0.6, 0.8), size = 3, 
                   straight = TRUE) +
  facet_wrap(~ Dominance, 
             scales = "free", nrow = 3, strip.position = c("top"),
             labeller = label_glue('({.l}) {`Dominance`}')) +
  labs(y = "_s~m~_", 
       x = "_s~f~_",
       fill = "_A~f~_ allele freq") +
  guides(fill = guide_colourbar(barheight = 15)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        text =  element_text(size = 12),
        axis.title.x = element_markdown(size = 18),
        axis.title.y = element_markdown(size = 18),
        legend.position = "right",
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6),
        strip.text.x = element_text(size = 11),
        legend.title = ggtext::element_markdown())

Figure S2: Predicted equilibrium frequency of the \(A_f\) allele, calculated from the population genetic model. The plot shows how varying the strength of the selection coefficient in females (\(s_f\)) and males (\(s_m\)) affects the fate of the \(A_f\) allele. The dashed lines enclose the parameter space where the \(A_f\) allele is predicted to stabilise at an intermediate frequency. Areas above and to the left of the upper dashed line indicate conditions where the allele is purged from the population, and areas below and to the right of the lower line indicate conditions where the allele spreads to fixation. Six dominance conditions are plotted. a additive expression, where \(h_f = h_m = 0.5\). b dominant expression of the \(A_f\) allele, where \(h_f = 1\) and \(h_m = 0\). c recessive expression of the \(A_f\) allele, where \(h_f = 0\) and \(h_m = 1\). d partial deleterious dominance reversal, where the \(A_f\) allele is recessive in females (\(h_f = 0.25\)) and dominant in males (\(h_m = 0.25\), where this is the dominance of the \(A_m\) allele). e and f show cases of beneficial dominance reversal, where the \(A_f\) allele is dominant when expressed in a female and recessive when expressed in a male. e shows an intermediate case where dominance is incomplete \(h_f = 0.75\) and \(h_m = 0.75\), while f depicts the parameter space for polymorphism when beneficial dominance reversal is complete (\(h_f = 1\) and \(h_m = 1\)).

As found previously, the conditions for polymorphism are very restrictive when selection is weak, but expand significantly when selection is strong and/or there is beneficial reversal of dominance i.e. when the \(A_f\) allele is (partially) dominant in females and (partially) recessive in males. Panel d shows why deleterious dominance reversal is not expected to be commonly observed in the genome. To contextualise \(s_f\) and \(s_m\), a value of 0.5 indicates that the focal allele is 50% fitter than its alternative in that particular environment, while a value of 1 indicates it has twice the fitness of the alternative.

\(~\)

Phase 2: modelling allele frequency change under sex-limited evolution

Calculating the per generation change in allele frequency

To predict the effect of sex-limited experimental evolution, we simulate the evolution of a female-beneficial, sexually antagonistic allele over 20 generations of female-limited inheritance. Our synthetic sex-limited inheritance regimes create haploid populations of autosomes, that are nearly always (see recombination compartment) expressed as heterozygotes with non-evolving, single genotype populations of balancer homologs, with which they cannot recombine. We therefore model the locus as haploid, with the exception that dominance relations can affect the strength of selection. The extent to which dominance matters is contingent upon whether the \(A_f\) allele is present at the locus on the non-evolving balancer. If so, the \(A_f\) allele will always be expressed, if not, dominance relations matter. We model the homozygous and heterozygous cases using two related equations.

To model evolution under female-limited inheritance when the evolving and non-evolving allele are homozygous, we use equation 1 presented in Dapper et al (2023). The \(A_f\) allele experiences conditions akin to a maternally inherited cytoplasmic allele. The change in allele frequency due to selection per generation is now

\[\begin{equation} \Delta p = \frac{s_fpq}{W_{f-limited}} \end{equation}\]

where \(W_{f-limited} = 1 + s_fp\).

We then model the effect of female-limited inheritance when the balancer does not carry the \(A_f\) allele. The new equation takes the form

\[\begin{equation} \Delta p = \frac{s_fh_fpq}{W_{f-limited}} \end{equation}\]

where \(W_{f-limited} = 1 + s_fh_fp\).

This equation helps point out that our design targets only those alleles that are expressed against the balancer genetic background. We therefore expect no response to selection for fully recessive alleles heterozygous with alternative alleles carried by the balancer chromosomes.

Modeling with or without dominance

Dominance doesn’t matter

Build the function

run_haploid_simulation_no_dominance <- function(gens_2, row, parameters){
  
  # Get the focal parameter space out of parameters[row], this makes it easier later
  
  generations <- gens_2
  prop_1 <- parameters$prop_1[row]
  s_f <- parameters$s_f[row]
  
  # Make the initial population
  pop <- 
    tibble(
      prop_2 = c(prop_1),
      W_female = 1 + prop_2*s_f)
  
  
  current_generation <- 1
  
  # Iterate over generations
  
  while(current_generation <= generations){
    
    
    pop <- pop %>%
      mutate(W_female = 1 + prop_2*s_f, # mean fitness changes as allele freqs change
             prop_2 = prop_2 + (s_f*prop_2*(1 - prop_2))/W_female, # here's the effect of selection
             generations_2 = gens_2)
    
    # Otherwise, increment generation counter by 1 and return to top of the 'while' loop
    
    current_generation <- current_generation + 1
  }
  
  pop %>% bind_cols(parameters[row, ]) # returns the final pop
  
}

Define the relevant parameter space

parameters_haploid_no_dominance <- 
  expand_grid(prop_1 = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9),
              s_f = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1))

Run the haploid simulation

haploid_results_no_dominance <- 
  run_haploid_simulation_no_dominance(gens_2 = 1, 1:nrow(parameters_haploid_no_dominance), parameters_haploid_no_dominance) %>% 
  mutate(generations_2 = 0, prop_2 = prop_1) %>% 
  bind_rows(map_dfr(1:20, run_haploid_simulation_no_dominance, 1:nrow(parameters_haploid_no_dominance), parameters_haploid_no_dominance)) %>%
  mutate(selection = as.factor(s_f),
         `Initial allele frequency` = as.factor(prop_1))

Dominance matters

Build the function

run_haploid_simulation_dominance <- function(gens_2, row, parameters){
  
  # Get the focal parameter space out of parameters[row], this makes it easier later
  
  generations <- gens_2
  prop_1 <- parameters$prop_1[row]
  s_f <- parameters$s_f[row]
  h_f <- parameters$h_f[row]
  
  
  # Make the initial population
  pop <- 
    tibble(
      prop_2 = c(prop_1),
      W_female = 1 + prop_2*s_f*h_f)
  
  
  current_generation <- 1
  
  # Iterate over generations
  
  while(current_generation <= generations){
    
    
    pop <- pop %>%
      mutate(W_female = 1 + prop_2*s_f*h_f, # mean fitness changes as allele freqs change
             prop_2 = prop_2 + (s_f*h_f*prop_2*(1 - prop_2))/W_female, # here's the effect of selection
             generations_2 = gens_2)
    
    # Otherwise, increment generation counter by 1 and return to top of the 'while' loop
    
    current_generation <- current_generation + 1
  }
  
  pop %>% bind_cols(parameters[row, ]) # returns the final pop
  
}

Define the relevant parameter space

parameters_haploid_dominance <- 
  expand_grid(prop_1 = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9),
              s_f = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),
              h_f = c(0, 0.25, 0.5, 0.75, 1))

Run the haploid simulation

haploid_results_dominance <- 
  run_haploid_simulation_dominance(gens_2 = 1, 1:nrow(parameters_haploid_dominance), parameters_haploid_dominance) %>% 
  mutate(generations_2 = 0, prop_2 = prop_1) %>% 
  bind_rows(map_dfr(1:20, run_haploid_simulation_dominance, 1:nrow(parameters_haploid_dominance), parameters_haploid_dominance)) %>%
  mutate(selection = as.factor(s_f),
         `Initial allele frequency` = as.factor(prop_1))

Plot the results

h_1 <-
  haploid_results_no_dominance %>%
  filter(`Initial allele frequency` != 0.9) %>% 
  ggplot(aes(x = generations_2, y = prop_2, group = selection, colour = s_f)) +
  geom_line(linewidth = 0.8) +
  scale_colour_gradientn(colours = pal1, limits = c(0.01, 1), 
                         breaks = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)) +
  guides(colour = guide_colourbar(barheight = 20)) +
  facet_wrap(~`Initial allele frequency`, labeller = label_glue(
    'Initial _A~f~_ allele freq: {`Initial allele frequency`}')) +
  labs(x = "Generations of female-limited inheritance",
       y = "Frequency of the female-benefiting (_A~f~_) allele",
       title = "(a) Allele always expressed",
       colour = "_s~f~_") +
  scale_x_continuous(breaks = c(0, 5, 10, 15, 20), limits = c(0, 20)) +
  theme_bw() +
  theme(axis.title.y = element_markdown(size = 14),
        axis.title.x = element_text(size = 14),
        legend.title = element_markdown(size = 14),
        strip.background = element_rect(fill = "Aliceblue", linewidth = .5),
        strip.text = element_markdown(),
        legend.margin = margin(6, 6, 6, 6),
        plot.title = element_markdown())

h_2 <-  
  haploid_results_dominance %>% 
  filter(h_f == 1,
         `Initial allele frequency` != 0.9) %>% 
  ggplot(aes(x = generations_2, y = prop_2, group = selection, colour = s_f)) +
  geom_line(linewidth = 0.8) +
  scale_colour_gradientn(colours = pal1, limits = c(0.01, 1), 
                         breaks = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)) +
  guides(colour = guide_colourbar(barheight = 20)) +
  facet_wrap(~`Initial allele frequency`, labeller = label_glue(
    'Initial _A~f~_ allele freq: {`Initial allele frequency`}')) +
  labs(x = "Generations of female-limited inheritance",
       y = "Frequency of the female-benefiting (_A~f~_) allele",
       title = "(b) _h~f~_ = 1",
       colour = "_s~f~_") +
  scale_x_continuous(breaks = c(0, 5, 10, 15, 20), limits = c(0, 20)) +
  theme_bw() +
  theme(axis.title.y = element_markdown(size = 14),
        axis.title.x = element_text(size = 14),
        legend.title = element_markdown(size = 14),
        strip.background = element_rect(fill = "Aliceblue", linewidth = .5),
        strip.text = element_markdown(),
        legend.margin = margin(6, 6, 6, 6),
        plot.title = element_markdown())

h_3 <-
  haploid_results_dominance %>% 
  filter(h_f == 0.5,
         `Initial allele frequency` != 0.9) %>% 
  ggplot(aes(x = generations_2, y = prop_2, group = selection, colour = s_f)) +
  geom_line(linewidth = 0.8) +
  scale_colour_gradientn(colours = pal1, limits = c(0.01, 1), 
                         breaks = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)) +
  guides(colour = guide_colourbar(barheight = 20)) +
  facet_wrap(~`Initial allele frequency`, labeller = label_glue(
    'Initial _A~f~_ allele freq: {`Initial allele frequency`}')) +
  labs(x = "Generations of female-limited inheritance",
       y = "Frequency of the female-benefiting (_A~f~_) allele",
       title = "(c ) _h~f~_ = 0.5",
       colour = "_s~f~_") +
  scale_x_continuous(breaks = c(0, 5, 10, 15, 20), limits = c(0, 20)) +
  theme_bw() +
  theme(axis.title.y = element_markdown(size = 14),
        axis.title.x = element_text(size = 14),
        legend.title = element_markdown(size = 14),
        strip.background = element_rect(fill = "Aliceblue", linewidth = .5),
        strip.text = element_markdown(),
        legend.margin = margin(6, 6, 6, 6),
        plot.title = element_markdown())

h_4 <-
  haploid_results_dominance %>% 
  filter(h_f == 0.25,
         `Initial allele frequency` != 0.9) %>% 
  ggplot(aes(x = generations_2, y = prop_2, group = selection, colour = s_f)) +
  geom_line(linewidth = 0.8) +
  scale_colour_gradientn(colours = pal1, limits = c(0.01, 1), 
                         breaks = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)) +
  guides(colour = guide_colourbar(barheight = 20)) +
  facet_wrap(~`Initial allele frequency`, labeller = label_glue(
    'Initial _A~f~_ allele freq: {`Initial allele frequency`}')) +
  labs(x = "Generations of female-limited inheritance",
       y = "Frequency of the female-benefiting (_A~f~_) allele",
       title = "(d) _h~f~_ = 0.25",
       colour = "_s~f~_") +
  scale_x_continuous(breaks = c(0, 5, 10, 15, 20), limits = c(0, 20)) +
  theme_bw() +
  theme(axis.title.y = element_markdown(size = 14),
        axis.title.x = element_text(size = 14),
        legend.title = element_markdown(size = 14),
        strip.background = element_rect(fill = "Aliceblue", linewidth = .5),
        strip.text = element_markdown(),
        legend.margin = margin(6, 6, 6, 6),
        plot.title = element_markdown())

h_1 + h_2 + h_3 + h_4 +
  plot_layout(guides = "collect")

Figure S3: The increase in frequency of the female-beneficial allele across 20 generations of female-limited inheritance. Curves show how varying the female selection coefficient (\(s_f\); increasing by increments of 0.1) for the allele changes the evolutionary trajectory, while comparing panels shows how this is affected by the initial \(A_f\) allele frequency. a the allele is always expressed, as would occur when the balancer homolog allele is also the \(A_f\) allele. b-d show cases where the balancer genotype carries the alternative allele to the \(A_f\) allele, where dominance relations matter.

Figure S3 shows that when selection is strong (\(s_f > 0.4\)), the \(A_f\) allele is very likely to reach fixation after 20 generations, except when the allele is rarely expressed. Even when \(s_f = 0.1\), the \(A_f\) allele can reach high frequencies after 20 generations of sex-limited evolution. The figure also illustrates that allele frequency change is non-linear. At intermediate frequencies change is fastest, whereas it is much slower at the floor or ceiling. Therefore, evolution from intermediate to high frequencies should result in greater change than evolution from high to intermediate frequencies.

\(~\)

Phase 3: can alleles return to their original equilibrium frequencies?

Now we can return to the phase 1 equations to assess the likelihood of alleles evolving back to their pre-experimental evolution frequencies. This time, we run the simulation for 20 generations and calculate \(\Delta A_f\), the change in allele frequency from generation 0 to 20 of relaxed biparental inheritance conditions that were enforced prior to the fitness assays.

Define a new parameter space

diploid_parameters_2 <- 
  expand_grid(
    s_f = c(0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1),
    s_m = c(seq(0, 1, length = resolution*4)),
    h_f = c(0, 0.25, 0.5, 0.75, 1),
    h_m = c(0, 0.25, 0.5, 0.75, 1),
    initial_freq_F =  c(seq(0.01, 1, length = resolution*2)))

Run the simulation

The populations experienced these conditions for 20 generations, which we can plug into the function.

diploid_rebound_results <- 
    run_diploid_simulation(20, 1:nrow(diploid_parameters_2), diploid_parameters_2)

Plot the results

plotting_data <-
  diploid_rebound_results %>% 
  mutate(F_change = prop_1 - initial_freq_F,
         Dominance = case_when(h_f == 0.5 & h_m == 0.5 ~ "Additive",
                               h_f == 1 & h_m == 0 ~ "Dominant",
                               h_f == 0 & h_m == 1 ~ "Recessive",
                               h_f == 0.75 & h_m == 0.75 ~ "Dominance reversal",
                               #h_f == 0.75 & h_m == 0.25 ~ "Deleterious dominance reversal",
                               .default = NA)) %>% 
  filter(!is.na(Dominance),
         s_f < 0.6) %>% 
  mutate(Dominance = fct_relevel(Dominance, 
                                 "Additive", "Dominant", "Recessive", 
                                 "Dominance reversal"))
  
  

#rebound_p1 <-
plotting_data %>%
  filter(s_f != 0.01) %>% 
  ggplot(aes(x = initial_freq_F, y = s_m)) +
  geom_blank() +
  geom_raster(aes(fill = F_change), alpha = 0.75) +
  scale_fill_gradientn(colours = pal1, breaks = c(-1, -0.5, 0, 0.5, 1), limits = c(-1, 1)) +
  #stat_contour(aes(z = F_change), colour = "black", linetype = 1,
             # breaks = c(-0.1, -0.2, -0.3, -0.4, -0.5, -0.6, -0.7, -0.8, -0.9, -0.99, 0.1, 0.2, 0.3, 0.4, 0.5)) +
  geom_textcontour(aes(z = F_change), 
                   breaks = c(-0.1, -0.2, -0.3, -0.4, -0.5, -0.6, -0.7, -0.8, -0.9, -0.99, 0.1, 0.2, 0.3, 0.4, 0.5),
                   size = 3, straight = TRUE) +
  stat_contour(aes(z = F_change), colour = "black", breaks = 0,
               linetype = 2, linewidth = 1) +
  
  facet_wrap(s_f ~ Dominance, 
             scales = "free", ncol = 4, strip.position = c("top"),
             labeller = label_glue('{`Dominance`}, _s~f~_ = {`s_f`}')) +
  labs(y = expression(s[m]), # this might seem backwards,
       x = "_A~f~_ freq after experimental evolution",
       fill = "Rebound change in _A~f~_ frequency") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) +
  guides(fill = guide_colourbar(barwidth = 10)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        text =  element_text(size = 12),
        axis.title = element_text(size = 18),
        axis.title.x = element_markdown(size = 18),
        legend.title = element_markdown(size = 18),
        legend.text = element_text(size = 12),
        legend.position = "bottom",
        legend.margin = margin(6, 6, 6, 6),
        strip.text.x = element_markdown(size = 13))

Figure S4: Predicted rebound in the equilibrium frequency of the \(A_f\) allele, after 20 generations of selection in both sexes, calculated from the population genetic model. The plot shows how varying the strength of the selection coefficient in males (\(s_m\)) and the frequency of the \(A_F\) allele at the end of experimental evolution affect the fate of the \(A_F\) allele. The dashed line indicates the conditions where no change is predicted for the frequency of the \(A_F\) allele. Contours show increments of 0.1 change in frequency, with warm regions indicating positive change and cool regions indicating a decline in frequency. Four dominance conditions are plotted across the columns: additivity, complete dominance, complete recessiveness and beneficial dominance reversal, where the \(A_F\) allele is partially dominant when expressed in a female (\(h_f = 0.75\)) and partially recessive when expressed in a male (\(h_m = 0.75\), where this refers to the dominance of the M allele in males). The selection coefficient in females (\(s_f\)) increases down the rows. Only cases where \(s_f\) tangibly changes allele frequencies while plausibly maintaining polymorphism during the sex-limited evolution phase are plotted.

Figure S4 shows that the \(A_f\) allele can decline in frequency when \(s_m\) is large relative to \(s_f\). However, as indicated by the dashed line’s position on the x axis in many of the panels, when selection is relatively equal between the sexes, very little change in frequency is predicted. Even when the \(A_m\) allele is dominant in males (e.g. Figure S4 columns 3-4), movement from high to intermediate frequencies requires \(s_m\) to be substantially stronger than \(s_m\).

\(~\)

Integrate the phases into a single simulation

Tweak the diploid simulation function

We need to tweak the run_diploid_simulation function to take the appropriate parameter values calculated in the prior simulation phases. The new run_diploid_simulation_2 function is identical, except it now accepts parameters with new names e.g. gens_3 instead of gens_1 and prop_2 instead of initial_freq_F. It also generates parameters with new names e.g. generations_3 and prop_3.

run_diploid_simulation_2 <- function(gens_3, row, parameters){
  
  
  generations_3 <- gens_3
  s_f <- parameters$s_f[row]
  s_m <- parameters$s_m[row]
  h_f <- parameters$h_f[row]
  h_m <- parameters$h_m[row]
  prop_2 <- parameters$prop_2[row]
  
  
  # Make the initial population
  pop <- 
    tibble(prop_3 = c(prop_2),
           generations_3 = gens_3)
  
  
  current_generation <- 1
  
  # Iterate over generations
  
  while(current_generation <= generations_3){
    pop <- pop %>%
     mutate(# here's the effect of selection
      f_freq_after_selection = ((prop_3^2)*(1 + s_f) + 0.5*(2*prop_3*(1 - prop_3))*(1 + s_f*h_f))/(
       (prop_3^2)*(1 + s_f) + (2*prop_3*(1 - prop_3))*(1 + s_f*h_f) + (1-prop_3)^2),
      m_freq_after_selection = ((prop_3^2) + 0.5*(2*prop_3*(1 - prop_3))*(1 + s_m*h_m))/(
       (prop_3^2) + (2*prop_3*(1 - prop_3))*(1 + s_m*h_m) + ((1-prop_3)^2)*(1 + s_m)),
      prop_3 = (f_freq_after_selection + m_freq_after_selection)/2,
     generations_3 = gens_3)
    
    # Otherwise, increment generation counter by 1 and return to top of the 'while' loop
    
    current_generation <- current_generation + 1
  }
  
  pop %>% bind_cols(parameters[row, ]) # returns the final pop
  
}

Build the combined function

combined_function <- function(parameters, gens_1, gens_2, gens_3){
  # phase 1
  
  polymorphism_results <- run_diploid_simulation(gens_1, 1:nrow(parameters), parameters) %>% 
    # just get the relevant cases we explore
    mutate(Dominance = case_when(h_f == 0.5 & h_m == 0.5 ~ "Additive",
                                 h_f == 0 & h_m == 1 ~ "Recessive (complete)",
                                 h_f == 0.25 & h_m == 0.75 ~ "Recessive (partial)",
                                 h_f == 1 & h_m == 0 ~ "Dominant",
                                 h_f == 0.25 & h_m == 0.25 ~ 
                                   "Deleterious dominance reversal (partial)",
                                 h_f == 0.75 & h_m == 0.75 ~ 
                                   "Beneficial dominance reversal (partial)",
                                 h_f == 1 & h_m == 1 ~ 
                                   "Beneficial dominance reversal (complete)",
                                 .default = NA)) %>% 
    filter(!is.na(Dominance)) %>% 
    mutate(prop_1 = case_when(prop_1 > 0.999 ~ 1,
                              prop_1 < 0.001 ~ 0,
                              .default = prop_1)) %>% 
    select(!contains("after")) # we don't need these to run the subsequent phases
  
  # phase 2
  
  haploid_results <- 
    run_haploid_simulation_dominance(gens_2, 1:nrow(polymorphism_results), 
                                     polymorphism_results) %>% 
    mutate(Dominance_matters = "YES") %>% 
    select(-W_female) %>% 
    bind_rows(run_haploid_simulation_no_dominance(gens_2, 1:nrow(polymorphism_results), 
                                                  polymorphism_results) %>%
                mutate(Dominance_matters = "NO") %>% 
                select(-W_female))
  
  # phase 3
  
  diploid_rebound_results <- 
    run_diploid_simulation_2(gens_3, 1:nrow(haploid_results), haploid_results)
  
  # order relevant columns and remove unnecessary ones for a nice output
  diploid_rebound_results %>% 
    select(-c(f_freq_after_selection , m_freq_after_selection )) %>% 
    select(initial_freq_F, s_f, s_m, h_f, h_m, Dominance, Dominance_matters, prop_1, prop_2, prop_3,
           generations_1, generations_2, generations_3, parameter_space_ID)
}

Run the combined function

To run the combined_function all we have to do is feed it the initial phase 1 parameter space and the generations we would like to run each phase for.

if(!file.exists("Integrated_results.csv")){
  Integrated_results <- combined_function(equilibrium_parameters,
                                          gens_1 = 10000, # phase 1
                                          gens_2 = 20, # phase 2
                                          gens_3 = 20) # phase 3
  write_csv(Integrated_results, file = "Integrated_results.csv")
} else Integrated_results <- read_csv("Integrated_results.csv")

Plot the results

I_1 <- 
  Integrated_results %>% 
  filter(Dominance_matters == "YES",
         Dominance != "Deleterious dominance reversal (partial)") %>% 
  mutate(Dominance = fct_relevel(Dominance, 
                                 "Additive", "Dominant",
                                 "Recessive (complete)",
                                 "Recessive (partial)",
                                 "Beneficial dominance reversal (partial)",
                                 "Beneficial dominance reversal (complete)")) %>% 
  
  ggplot(aes(x = s_f, y = s_m)) +
  geom_blank() +
  geom_raster(aes(fill = prop_3), alpha = 0.75) +
  scale_fill_gradientn(colours = pal1, breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)) +
  #stat_contour(aes(z = prop_3), colour = "black", linetype = 1,
   #           breaks = c(0.1, 0.25, 0.5, 0.75, 0.9)) +
  stat_contour(aes(z = prop_3), colour = "black", breaks = c(0.01, 0.99),
               linetype = 2, linewidth = 1) +
  
  facet_wrap( ~ Dominance, 
             scales = "free", ncol = 2, strip.position = c("top"),
             labeller = label_glue('{`Dominance`}')) +
  labs(y = "_s~m~_", 
       x = "_s~f~_",
       fill = "Final _A~f~_ allele frequency",
       title = "(a) Dominance matters in treatment phase") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) +
  guides(fill = guide_colourbar(barwidth = 10)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        text =  element_text(size = 11),
        axis.title.x = element_markdown(size = 20),
        axis.title.y = element_markdown(size = 20),
        legend.title = element_markdown(size = 14),
        legend.text = element_text(size = 12),
        legend.position = "bottom",
        legend.margin = margin(6, 6, 6, 6),
        strip.text.x = element_text(size = 11))

I_2 <-
  Integrated_results %>% 
  filter(Dominance_matters == "NO",
         Dominance != "Deleterious dominance reversal (partial)") %>% 
  mutate(Dominance = fct_relevel(Dominance, 
                                 "Additive", "Dominant",
                                 "Recessive (complete)",
                                 "Recessive (partial)",
                                 "Beneficial dominance reversal (partial)",
                                 "Beneficial dominance reversal (complete)")) %>% 
  
  ggplot(aes(x = s_f, y = s_m)) +
  geom_blank() +
  geom_raster(aes(fill = prop_3), alpha = 0.75) +
  scale_fill_gradientn(colours = pal1, breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)) +
  #stat_contour(aes(z = prop_3), colour = "black", linetype = 1,
  #           breaks = c(0.1, 0.25, 0.5, 0.75, 0.9)) +
  stat_contour(aes(z = prop_3), colour = "black", breaks = c(0.01, 0.99),
               linetype = 2, linewidth = 1) +
  
  facet_wrap( ~ Dominance, 
              scales = "free", ncol = 2, strip.position = c("top"),
              labeller = label_glue('{`Dominance`}')) +
  labs(y = "_s~m~_", 
       x = "_s~f~_",
       fill = "Final _A~f~_ allele frequency",
       title = "(b) Allele always expressed in treatment phase") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) +
  guides(fill = guide_colourbar(barwidth = 10)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        text =  element_text(size = 11),
        axis.title.x = element_markdown(size = 20),
        axis.title.y = element_markdown(size = 20),
        legend.title = element_markdown(size = 14),
        legend.text = element_text(size = 12),
        legend.position = "bottom",
        legend.margin = margin(6, 6, 6, 6),
        strip.text.x = element_text(size = 11))

I_1 + I_2 +
  plot_layout(guides = "collect") & theme(legend.position = 'bottom')

Figure S5: the final frequency of the \(A_f\) allele following sex-limited evolution and subsequent relaxation of sex-limited inheritance. The dashed lines bound the parameter space where polymorphism is still predicted. Note that compared with Figure S2, the parameter space for polymorphism is reduced, as the \(A_f\) increases to fixation in many cases during the sex-limited evolution phase. a the allele is always expressed during the sex-limited inheritance phase, as would occur when the balancer homolog allele is also the \(A_f\) allele. b-d show cases where the balancer genotype carries the alternative allele to the \(A_f\) allele, where dominance relations matter during the sex-limited inheritance phase.

I_3 <-
  Integrated_results %>% 
  mutate(Haploid_change = prop_2 - prop_1,
         Relaxed_change = prop_3 - prop_2,
         Total_change = prop_3 - prop_1) %>% 
  filter(Dominance_matters == "NO",
         Dominance != "Deleterious dominance reversal (partial)") %>% 
  mutate(Dominance = fct_relevel(Dominance, 
                                 "Additive", "Dominant",
                                 "Recessive (complete)",
                                 "Recessive (partial)",
                                 "Beneficial dominance reversal (partial)",
                                 "Beneficial dominance reversal (complete)")) %>% 
  ggplot(aes(x = s_f, y = s_m)) +
  geom_blank() +
  geom_raster(aes(fill = Total_change), alpha = 0.75) +
  scale_fill_gradientn(colours = pal1, breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1),
                       limits = c(-0.05, 1)) +
  stat_contour(aes(z = Total_change), colour = "black", linetype = 1,
               breaks = c(0.2, 0.4, 0.6, 0.8)) +
  #geom_textcontour(aes(z = Total_change), 
  #               breaks = c(0.2, 0.4, 0.6, 0.8),
  #              size = 2, straight = TRUE) +
  stat_contour(aes(z = Total_change), colour = "black", breaks = c(0.01),
               linetype = 2, linewidth = 1) +
  
  facet_wrap( ~ Dominance, 
              scales = "free", ncol = 2, strip.position = c("top"),
              labeller = label_glue('{`Dominance`}')) +
  labs(y = "_s~m~_", 
       x = "_s~f~_",
       fill = "Effect of experiment on _A~f~_ freq",
       title = "(a) Allele always expressed in treatment phase") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) +
  guides(fill = guide_colourbar(barwidth = 10)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        text =  element_text(size = 11),
        axis.title.x = element_markdown(size = 20),
        axis.title.y = element_markdown(size = 20),
        legend.title = element_markdown(size = 14),
        legend.text = element_text(size = 12),
        legend.position = "bottom",
        legend.margin = margin(6, 6, 6, 6),
        strip.text.x = element_text(size = 11))



I_4 <-
  Integrated_results %>% 
  mutate(Haploid_change = prop_2 - prop_1,
         Relaxed_change = prop_3 - prop_2,
         Total_change = prop_3 - prop_1) %>% 
  filter(Dominance_matters == "YES",
         Dominance != "Deleterious dominance reversal (partial)") %>% 
  mutate(Dominance = fct_relevel(Dominance, 
                                 "Additive", "Dominant",
                                 "Recessive (complete)",
                                 "Recessive (partial)",
                                 "Beneficial dominance reversal (partial)",
                                 "Beneficial dominance reversal (complete)")) %>% 
  ggplot(aes(x = s_f, y = s_m)) +
  geom_blank() +
  geom_raster(aes(fill = Total_change), alpha = 0.75) +
  scale_fill_gradientn(colours = pal1, breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1),
                       limits = c(-0.05, 1)) +
  stat_contour(aes(z = Total_change), colour = "black", linetype = 1,
               breaks = c(0.2, 0.4, 0.6, 0.8)) +
  #geom_textcontour(aes(z = Total_change), 
  #               breaks = c(0.2, 0.4, 0.6, 0.8),
  #              size = 2, straight = TRUE) +
  stat_contour(aes(z = Total_change), colour = "black", breaks = c(0.01),
               linetype = 2, linewidth = 1) +
  facet_wrap( ~ Dominance, 
              scales = "free", ncol = 2, strip.position = c("top"),
              labeller = label_glue('{`Dominance`}')) +
  labs(y = "_s~m~_", 
       x = "_s~f~_",
       fill = "Effect of experiment on _A~f~_ freq",
       title = "(b) Dominance matters in treatment phase") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) +
  guides(fill = guide_colourbar(barwidth = 10)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        text =  element_text(size = 11),
        axis.title.x = element_markdown(size = 20),
        axis.title.y = element_markdown(size = 20),
        legend.title = element_markdown(size = 14),
        legend.text = element_text(size = 12),
        legend.position = "bottom",
        legend.margin = margin(6, 6, 6, 6),
        strip.text.x = element_text(size = 11))



I_3 + I_4 +
  plot_layout(guides = "collect") & theme(legend.position = 'bottom')

Figure S6: The effects of sex-limited inheritance for 20 generations, followed by 16 generations where a selection response could occur in either sex. The conditions where the \(A_f\) allele is expected to increase in frequency are bounded by the dashed lines, outside of which change = 0. Contour lines enclosed within the dashed lines indicate intervals of 0.2, increasing with the warmth of the colour displayed. a the allele is always expressed during the sex-limited inheritance phase, as would occur when the balancer homolog allele is also the \(A_f\) allele. b-d show cases where the balancer genotype carries the alternative allele to the \(A_f\) allele, where dominance relations matter during the sex-limited inheritance phase.

The combined simulation predicts that the relaxed selection phase should allow a partial shift back to previous frequencies for a subset of sexually antagonistic loci, but that the majority should retain a signature of sex-limited evolution. Areas of no change occur because the allele is predicted to either be fixed or purged from the population prior to the start of the experiment.

\(~\)

References

Kidwell, J. F., M. T. Clegg, F. M. Stewart, and T. Prout. 1977. Regions of stable equilibria for models of differential selection in the two sexes under random mating. Genetics 85:171–183.

Rice, W. R. 1984. Sex Chromosomes and the Evolution of Sexual Dimorphism. Evolution 38:735–742.

Curtsinger, J. W., Service, Philip M., and T. Prout. 1994. Antagonistic Pleiotropy, Reversal of Dominance, and Genetic Polymorphism. American Naturalist 144:210–228.

Connallon, T., R. M. Cox, and R. Calsbeek. 2009. Fitness consequences of sex-specific selection. Evolution 64:1671–1682.

Patten, M. M., D. Haig, and F. Ubeda. 2010. Fitness variation due to sexual antagonism and linkage disequilibrium. Evolution 64:3638–3642.

Connallon, T., and S. F. Chenoweth. 2019. Dominance reversals and the maintenance of genetic variation for fitness. PLoS Biology 17:e3000118.

Dapper, A. L., A. E. Diegel, and M. J. Wade. 2023. Relative rates of evolution of male-beneficial nuclear compensatory mutations and male-harming Mother’s Curse mitochondrial alleles. Evolution 77:1945-1955

\(~\)

Analysis

Load in the data

\(~\)

fitness_data <- read_csv("data/SLC_fitness_data.csv") %>% 
  mutate(Fitness_vial_ID = as.factor(Fitness_vial_ID),
         Block = as.factor(Block),
         Population = as.factor(Population),
         Treatment = as.factor(Treatment),
         GFP = as.factor(GFP),
         Sex = as.factor(Sex),
         Rearing_vial = as.factor(Rearing_vial),
         Total_red_offspring = Red_female_offspring + Red_male_offspring,
         Total_bw_offspring = bw_female_offspring + bw_male_offspring,
         Total_offspring = Total_red_offspring + Total_bw_offspring) %>% 
  rename(Inheritance_treatment = Treatment)

# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'fitness_data')),
      pageLength = 692
    )
  )
}


my_data_table(fitness_data %>% select(-Comment))

Column explanations

Fitness_vial_ID: a unique identifier for each trial of the fitness assay.

Block: the experiment was run in three distinct blocks, using flies from separate generations.

Population: we measured the fitness of flies from 12 independent populations that contained autosomes that had undergone experimental evolution.

Inheritance_treatment: the populations carried autosomes that had been exposed to one of three inheritance treatments for 20 generations: a female-limited inheritance treatment where autosomes were always passed from mother to daughter, a male-limited treatment where autosomes were passed from father to son, and a control condition where inheritance was unconstrained.

GFP: the GFP marker carried by the population. UBI indicates the presence of a transgene that encodes ubiquitous expression of GFP, while 3xP indicates the presence of a different transgene that encodes the expression of GFP in the ocelli.

Sex: the sex of the individuals that we were measuring the fitness of.

Rearing_vial: the vial the treatment flies used in the trial developed in. This variable is included to capture variation explained by the rearing environment e.g. small differences in food moisture content or quantity. Note that females and males can have the same rearing vial as the sexes were reared together.

Red_female_offspring: the number of adult female offspring sired/produced by flies sourced from one of the 12 populations.

Red_male_offspring: the number of adult male offspring sired/produced by flies sourced from one of the 12 populations.

bw_female_offspring: the number of adult female offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

bw_male_offspring: the number of adult male offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

Total_red_offspring: the total number (sexes pooled) of adult offspring sired/produced by flies sourced from one of the 12 populations.

Total_bw_offspring: the total number (sexes pooled) of adult offspring sired/produced by competitor flies.

Total_offspring: the total number (sexes and eye colours pooled) of adult offspring counted in each vial.

\(~\)

Modelling approach

\(~\)

Female and male fitness are fundamentally different concepts / traits. There are also several differences between our female and male fitness assays. The major difference is that the male assay contains half the number of females in any given vial than does the female assay. The logic behind this design choice is that sexually selected processes are a more important determinant of male fitness than they are female fitness, so any fitness differences may only be observed when competition for fertilisations is high.

For these reasons, we choose to split the data up and model female and male fitness separately.

female_fitness <- 
  fitness_data %>%
  filter(Sex == "Female")

male_fitness <- 
  fitness_data %>%
  filter(Sex == "Male") %>% 
  mutate(prop_red = Total_red_offspring / Total_offspring)

We fit the following fixed and random effects to model female and male fitness. Our aim is to estimate the causal effect that Inheritance_treatment (I) has on fitness (F).

Fixed effects

Inheritance_treatment (I): this is the inheritance regime that the autosomes carried by each of the populations were subject to for 20 generations. There are three levels: populations carrying female-adapted autosomes, populations containing male-adapted autosomes and populations carrying control autosomes that experienced an unmanipulated inheritance regime. We are designing our model to test for a causal effect of this variable.

Mediator variables

Block (B): fitness might differ between the three distinct blocks we split our experiment up into. Blocks differed temporally, used flies from different generations and different batches of food. It is also possible that there were minor fluctuations in the lighting and temperature environment experienced during development between blocks. Each of these variables may introduce variation into our fitness measurements, that can be accounted for by including the Block variable in our model.

GFP (G): it is possible that fitness may be affected by the GFP transgene carried by each population. For example, one could imagine that any unintended fitness effects of a transgene might be of greater magnitude if it is expressed in a larger proportion of tissues, as is the case for the UBI transgene versus the 3xP transgene. Note that each GFP type is carried by an equal number of populations from each of the three evolutionary treatments.

Varying/Random effects

Population (P): our design contained 12 independent populations of autosomes that originated from a single outbred laboratory fly population. The populations were split and autosomes from each were subjected to one of the three evolution treatments for 20 generations. 4 populations experienced a female-limited inheritance regime, 4 a male-limited regime and 4 an unlimited or control regime.

Rearing_vial (R): the vial individual flies developed within may introduce further variation into our response variable. Like Block this variable controls for micro-environmental variation.

\(~\)

Accounting for over-dispersion

The data is over-dispersed with several highly influential (outlier) observations that have large effects on our posterior prediction. To combat this, we fit models following the betabinomial distribution family, as this is better equipped to deal with extreme values at the tails i.e. overdispersion.

However, the beta-binomial is not a native family in brms, we need to create the distribution family using the custom_family() function. The code below is taken directly from the brms_customfamilies vignette, which can be viewed here.

beta_binomial2 <- custom_family(
  "beta_binomial2", dpars = c("mu", "phi"),
  links = c("logit", "log"), lb = c(NA, 2), 
  # note that we set the lower bound to 2, following McElreath, rather than Buerkner. This means that the most conservative estimate for phi we get is a flat expectation between 0 and 1
  type = "int", vars = "vint1[n]"
)


stan_funs <- "
  real beta_binomial2_lpmf(int y, real mu, real phi, int T) {
    return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);
  }
  int beta_binomial2_rng(real mu, real phi, int T) {
    return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);
  }
"

stanvars <- stanvar(scode = stan_funs, block = "functions")

\(~\)

Female fitness

\(~\)

To estimate the fitness of females carrying each of the three autosome types, we placed three experimental females into a yeasted vial with three female competitors that carried the bw mutation. We then introduced six males that also carried the bw mutation. We allowed them to mate and oviposit for three days, then removed all adults from the vial. 12 days later we counted all of the adult progeny in the vial and scored them for eye-colour. Progeny with red eyes were produced by the experimental females ( bw is recessive) and progeny with brown eyes were produced by the competitor females. We calculated fitness as the proportion of red eyed offspring in the vial.

We present the fixed effects from the model:

\(~\)

female_fitness_model <-
  brm(Total_red_offspring | vint(Total_offspring) ~ 1 + Inheritance_treatment + Block + GFP + (1|Population) + (1|Rearing_vial),
      data = female_fitness, family = beta_binomial2, 
      prior = c(prior(normal(0, 1), class = Intercept), # this sets the prior mean at 0.5 relative fitness
                prior(normal(0, 1.5), class = b),
                prior(exponential(1), class = phi)),
      iter = 8000, warmup = 4000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.98, max_treedepth = 10),
      seed = 2, stanvars = stanvars, file = "Fits/female_fitness.model")

fixef(female_fitness_model) %>% 
  kable(digits = 3) %>% 
  kable_styling()

	Estimate	Est.Error	Q2.5	Q97.5
Intercept	1.110	0.095	0.921	1.296
Inheritance_treatmentFemale_limited	0.174	0.104	-0.031	0.383
Inheritance_treatmentMale_limited	0.187	0.104	-0.014	0.398
Block2	-0.637	0.079	-0.792	-0.480
Block3	-0.624	0.071	-0.765	-0.487
GFPUBI	-0.323	0.084	-0.488	-0.156

We need to write some additional code to get some post processing stuff i.e. LOO to work. Code courtesy of the brms_customfamilies vignette, which can be viewed here.

Run LOO to see if we’ve effectively modelled over dispersion.

female_fitness_model <- add_criterion(female_fitness_model, criterion = "loo", file = "Fits/female_fitness.model")

loo(female_fitness_model)

## 
## Computed from 16000 by 327 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1400.2 10.7
## p_loo        23.1  2.1
## looic      2800.4 21.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

The beta-binomial model looks good. It returns no points with high pareto k values.

Conduct a posterior predictive check to confirm our model is doing what we want it to.

pp_check(female_fitness_model, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

The posterior predictive distribution recapitulates our raw data quite well, indicating the model is making sensible predictions.

\(~\)

Derive predictions from the posterior

Get predictions from the model

Table S2. Estimated female fitness for flies carrying autosomes derived from each of the three inheritance regimes.

draws <- as_draws_df(female_fitness_model)

draws_female <-
  draws  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(`Female-limited`, Control, `Male-limited`) %>% 
    pivot_longer(cols = c(`Female-limited`, Control, `Male-limited`),
                 names_to = "Inheritance treatment") %>% 
  rename(prop_focal_offspring = value) %>% 
  arrange(`Inheritance treatment`)

draws_female %>% 
  group_by(`Inheritance treatment`) %>% 
  summarise(`Estimated prop. of offspring produced` = mean(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance treatment	Estimated prop. of offspring produced	2.5%	97.5%
Control	0.752	0.715	0.785
Female-limited	0.783	0.749	0.814
Male-limited	0.785	0.752	0.815

Table S3. Differences in female fitness between each of the three inheritance regimes.

# Find differences and visualise in table

draws_diff <- draws  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>%
  mutate(`Female - Control` = `Female-limited` - Control,
         `Male - Control` = `Male-limited` - Control,
         `Female - Male` = `Female-limited` - `Male-limited`) %>% 
  select(`Female - Control`, `Male - Control`, `Female - Male`)

rbind(
  draws_diff %>% 
    summarise(`Diff in offspring produced per 100` = mean(`Male - Control`)*100,
              `2.5%` = quantile(`Male - Control`, probs = 0.025) *100,
              `97.5%` = quantile(`Male - Control`, probs = 0.975) *100) %>% 
    mutate(Contrast = "Male inherited - Control"),
  
  draws_diff %>% 
    summarise(`Diff in offspring produced per 100` = mean(`Female - Control`)*100,
              `2.5%` = quantile(`Female - Control`, probs = 0.025) *100,
              `97.5%` = quantile(`Female - Control`, probs = 0.975) *100) %>% 
    mutate(Contrast = "Female inherited - Control"),
  
  draws_diff %>% 
    summarise(`Diff in offspring produced per 100` = mean(`Female - Male`)*100,
              `2.5%` = quantile(`Female - Male`, probs = 0.025) *100,
              `97.5%` = quantile(`Female - Male`, probs = 0.975) *100) %>% 
    mutate(Contrast = "Female inherited - Male inherited")
) %>% 
  select(Contrast, `Diff in offspring produced per 100`, `2.5%`, `97.5%`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 2)

Contrast	Diff in offspring produced per 100	2.5%	97.5%
Male inherited - Control	3.32	-0.25	7.11
Female inherited - Control	3.09	-0.56	6.84
Female inherited - Male inherited	-0.23	-3.7	3.2

The effects of the block and GFP predictors

Table S4. The effects of the fixed predictor variables on female fitness that are not directly related to intralocus sexual conflict. Female fitness measured in Block 1 was higher than that measured in Blocks 2 and 3. Females carrying autosomes marked with 3xP GFP had higher fitness than those expressing UBI GFP.

draws_diff_other <- draws  %>% 
  select(b_Intercept, b_Block2, b_Block3, b_GFPUBI) %>% 
  mutate(`Block 1, 3xP` = inv_logit_scaled(b_Intercept),
         `Block 2` = inv_logit_scaled(b_Intercept + b_Block2),
         `Block 3` = inv_logit_scaled(b_Intercept + b_Block3),
         UBI = inv_logit_scaled(b_Intercept + b_GFPUBI)) %>%
  mutate(`Block 1 - Block 2` = `Block 1, 3xP` - `Block 2`,
         `Block 1 - Block 3` = `Block 1, 3xP` - `Block 3`,
         `3xP - UBI` = `Block 1, 3xP` - UBI) %>% 
  select(`Block 1 - Block 2`, `Block 1 - Block 3`, `3xP - UBI`)

# Find differences and visualise in table

rbind(
draws_diff_other %>% 
  summarise(`Diff in offspring produced per 100` = mean(`Block 1 - Block 2`)*100,
            `2.5%` = quantile(`Block 1 - Block 2`, probs = 0.025) *100,
            `97.5%` = quantile(`Block 1 - Block 2`, probs = 0.975) *100) %>% 
  mutate(Contrast = "Block 1 - Block 2"),


draws_diff_other %>% 
  summarise(`Diff in offspring produced per 100` = mean(`Block 1 - Block 3`)*100,
            `2.5%` = quantile(`Block 1 - Block 3`, probs = 0.025) *100,
            `97.5%` = quantile(`Block 1 - Block 3`, probs = 0.975) *100) %>% 
  mutate(Contrast = "Block 1 - Block 3"),


draws_diff_other %>% 
  summarise(`Diff in offspring produced per 100` = mean(`3xP - UBI`)*100,
            `2.5%` = quantile(`3xP - UBI`, probs = 0.025) *100,
            `97.5%` = quantile(`3xP - UBI`, probs = 0.975) *100) %>% 
  mutate(Contrast = "3xP - UBI")
) %>% 
  select(Contrast, `Diff in offspring produced per 100`, `2.5%`, `97.5%`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 2)

Contrast	Diff in offspring produced per 100	2.5%	97.5%
Block 1 - Block 2	13.59	10.1	17.08
Block 1 - Block 3	13.27	10.35	16.32
3xP - UBI	6.48	3.11	9.84

\(~\)

Male fitness

\(~\)

To estimate the fitness of males carrying each of the three chromosome types, we conducted an experiment very similar to the female fitness assay. However, because male fitness has stronger covariance with fertilisation events than does female fitness, we conducted the male fitness assay with a 1:2 sex ratio (female:male) rather than the 1:1 ratio used in the female assay. This increases the strength of sexual selection and is a more appropriate way to expose differences in fitness between groups of males. As with the females, we calculated fitness as the proportion of red-eyed offspring in the vial.

male_fitness_model <-
  brm(Total_red_offspring | vint(Total_offspring) ~ 1 + Inheritance_treatment + Block + GFP + (1|Population) + (1|Rearing_vial),
      data = male_fitness, family = beta_binomial2, 
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, 1.5), class = b),
                prior(exponential(1), class = phi)),
      iter = 8000, warmup = 4000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.95, max_treedepth = 10),
      seed = 2, stanvars = stanvars, file = "Fits/male_fitness.model")


fixef(male_fitness_model) %>% 
  kable(digits = 3) %>% 
  kable_styling()

	Estimate	Est.Error	Q2.5	Q97.5
Intercept	0.761	0.233	0.291	1.217
Inheritance_treatmentFemale_limited	0.422	0.270	-0.126	0.970
Inheritance_treatmentMale_limited	0.107	0.265	-0.428	0.632
Block2	0.970	0.148	0.682	1.261
Block3	-0.035	0.141	-0.311	0.240
GFPUBI	-0.336	0.218	-0.771	0.103

Run LOO to see if we’ve effectively modelled over dispersion.

male_fitness_model <- add_criterion(male_fitness_model, criterion = "loo", file = "Fits/male_fitness.model")

loo(male_fitness_model)

## 
## Computed from 16000 by 360 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1561.9 21.4
## p_loo        23.1  1.9
## looic      3123.8 42.8
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     357   99.2%   5284      
##  (0.5, 0.7]   (ok)         2    0.6%   5066      
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   1    0.3%   8000      
## See help('pareto-k-diagnostic') for details.

There is one point having a large effect on the posterior. Upon inspection, this data point is not an unreasonable one and there is no cause to remove it from the dataset. It also does not change the causal effect of inheritance treatment on male fitness.

Conduct the posterior predictive check…

pp_check(male_fitness_model, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

Derive estimates from posterior

Get predictions from the model

Table S5. Estimated male fitness for flies carrying autosomes derived from each of the three inheritance regimes.

draws_m <- as_draws_df(male_fitness_model)

# predictions averaged over mediator variables

draws_male <-
  draws_m  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(Control, `Female-limited`, `Male-limited`) %>% 
    pivot_longer(cols = c(Control, `Female-limited`, `Male-limited`),
                 names_to = "Inheritance treatment") %>% 
  rename(prop_focal_offspring = value)


draws_male %>% 
  group_by(`Inheritance treatment`) %>% 
  summarise(`Estimated prop. of offspring sired` = mean(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance treatment	Estimated prop. of offspring sired	2.5%	97.5%
Control	0.679	0.572	0.771
Female-limited	0.763	0.67	0.84
Male-limited	0.702	0.598	0.792

Table S6. Differences in male fitness between each of the three inheritance regimes.

draws_diff_m <- draws_m %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>%
  mutate(`Female - Control` = `Female-limited` - Control,
         `Male - Control` = `Male-limited` - Control,
         `Female - Male` = `Female-limited` - `Male-limited`) %>% 
  select(`Female - Control`, `Male - Control`, `Female - Male`)

# Find differences and visualise in table

rbind(
  draws_diff_m %>% 
    summarise(`Diff in offspring produced per 100` = mean(`Female - Control`)*100,
              `2.5%` = quantile(`Female - Control`, probs = 0.025) *100,
              `97.5%` = quantile(`Female - Control`, probs = 0.975) *100) %>% 
    mutate(Contrast = "Female inherited - Control"),
  
  draws_diff_m %>% 
    summarise(`Diff in offspring produced per 100` = mean(`Male - Control`)*100,
              `2.5%` = quantile(`Male - Control`, probs = 0.025) *100,
              `97.5%` = quantile(`Male - Control`, probs = 0.975) *100) %>% 
    mutate(Contrast = "Male inherited - Control"),
  
  
  draws_diff_m %>% 
    summarise(`Diff in offspring produced per 100` = mean(`Female - Male`)*100,
              `2.5%` = quantile(`Female - Male`, probs = 0.025) *100,
              `97.5%` = quantile(`Female - Male`, probs = 0.975) *100) %>% 
    mutate(Contrast = "Female inherited - Male inherited")
) %>% 
  select(Contrast, `Diff in offspring produced per 100`, `2.5%`, `97.5%`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 2)

Contrast	Diff in offspring produced per 100	2.5%	97.5%
Female inherited - Control	8.33	-2.44	19.22
Male inherited - Control	2.26	-9.22	13.44
Female inherited - Male inherited	6.08	-4.42	16.8

The effects of the block and GFP predictors

Table S7. The effects of the fixed predictor variables on male fitness that are not directly related to intralocus sexual conflict. Male fitness measured in Block 2 was higher than that measured in Blocks 1 and 3. There was no effect of GFP transgene on male fitness.

draws_diff_other_m <- draws_m  %>% 
  select(b_Intercept, b_Block2, b_Block3, b_GFPUBI) %>% 
  mutate(`Block 1, 3xP` = inv_logit_scaled(b_Intercept),
         `Block 2` = inv_logit_scaled(b_Intercept + b_Block2),
         `Block 3` = inv_logit_scaled(b_Intercept + b_Block3),
         UBI = inv_logit_scaled(b_Intercept + b_GFPUBI)) %>%
  mutate(`Block 1 - Block 2` = `Block 1, 3xP` - `Block 2`,
         `Block 1 - Block 3` = `Block 1, 3xP` - `Block 3`,
         `3xP - UBI` = `Block 1, 3xP` - UBI) %>% 
  select(`Block 1 - Block 2`, `Block 1 - Block 3`, `3xP - UBI`)

# Find differences and visualise in table

rbind(
draws_diff_other_m %>% 
  summarise(`Diff in offspring produced per 100` = mean(`Block 1 - Block 2`)*100,
            `2.5%` = quantile(`Block 1 - Block 2`, probs = 0.025) *100,
            `97.5%` = quantile(`Block 1 - Block 2`, probs = 0.975) *100) %>% 
  mutate(Contrast = "Block 1 - Block 2"),


draws_diff_other_m %>% 
  summarise(`Diff in offspring produced per 100` = mean(`Block 1 - Block 3`)*100,
            `2.5%` = quantile(`Block 1 - Block 3`, probs = 0.025) *100,
            `97.5%` = quantile(`Block 1 - Block 3`, probs = 0.975) *100) %>% 
  mutate(Contrast = "Block 1 - Block 3"),


draws_diff_other_m %>% 
  summarise(`Diff in offspring produced per 100` = mean(`3xP - UBI`)*100,
            `2.5%` = quantile(`3xP - UBI`, probs = 0.025) *100,
            `97.5%` = quantile(`3xP - UBI`, probs = 0.975) *100) %>% 
  mutate(Contrast = "3xP - UBI")
) %>% 
  select(Contrast, `Diff in offspring produced per 100`, `2.5%`, `97.5%`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 2)

Contrast	Diff in offspring produced per 100	2.5%	97.5%
Block 1 - Block 2	-16.76	-23.17	-10.91
Block 1 - Block 3	0.75	-5.25	6.8
3xP - UBI	7.62	-2.31	17.49

\(~\)

Building Figure 2

# female plots

f_1 <-
  draws_female %>% 
  mutate(`Inheritance treatment` = fct_relevel(`Inheritance treatment`, "Female-limited", "Control", "Male-limited")) %>% 
  ggplot(aes(`Inheritance treatment`, prop_focal_offspring)) + 
  stat_halfeye(fill = "grey", .width = c(0.66, 0.95), alpha = 1,
               point_interval = "mean_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5) + # width indicates the uncertainty intervals: we have 66% and 95% intervals
  #scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  coord_flip() +
  ylab("Female fitness\n(prop. offspring produced)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_2 <-
  draws_diff %>% 
  gather(key = parameter, value = `Fitness difference`) %>% 
  as_tibble() %>%   

  ggplot(aes(parameter, `Fitness difference`)) + 
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, point_interval = "mean_qi",
               slab_fill = "grey",
               shape = 21, point_size = 3, stroke = 1.5,
               point_fill = "white") + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab("Treatment contrast") +
  ylab("Female fitness difference\n(prop. offspring produced)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

# male plots

f_3 <-
  draws_male %>% 
  mutate(`Inheritance treatment` = fct_relevel(`Inheritance treatment`, "Female-limited", "Control", "Male-limited")) %>% 
  ggplot(aes(`Inheritance treatment`, prop_focal_offspring)) + 
   stat_halfeye(fill = "grey", .width = c(0.66, 0.95), alpha = 1,
               point_interval = "mean_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  coord_flip() +
  ylab("Male fitness\n(prop. offspring produced)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_4 <-
  draws_diff_m %>%
  gather(key = parameter, value = `Fitness difference`) %>% 
  as_tibble() %>% 
  ggplot(aes(parameter, `Fitness difference`)) + 
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, point_interval = "mean_qi",
               slab_fill = "grey",
               shape = 21, point_size = 3, stroke = 1.5,
               point_fill = "white") + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hokusai3", 3)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(breaks = c(-0.2, -0.1, 0, 0.1, 0.2)) +
  xlab("Treatment contrast") +
  ylab("Male fitness difference\n(prop. offspring sired)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())


(f_1 + f_2) /(f_3 + f_4) + 
  plot_annotation(tag_levels = 'a')

Figure 2: a shows the estimated distribution of the mean for female fitness for flies carrying autosomes that had previously experienced unconstrained inheritance (control), female-limited inheritance or male-limited inheritance. b shows the difference contrast in female fitness between each of the three inheritance treatments. This difference is on the proportion scale, where a value of 0.1 indicates that females of a given inheritance treatment produce 10 more offspring per every 100 when cohabiting with bw competitor females. c and d depict the same things as a and b, except for male fitness.